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We describe BlackHat, an automated C++ program for calculating one-loop amplitudes, and the techniques 
used in its construction. These include the unitarity method and on-shell recursion. The other ingredients 
are compact analytic formulae for tree amplitudes for four- dimensional helicity states. The program computes 
amplitudes numerically, using analytic formula only for the tree amplitudes, the starting point for the recursion, 
and the loop integrals. We make use of recently developed on-shell methods for evaluating coefficients of loop 
integrals, in particular a discrete Fourier projection as a means of improving numerical stability. We illustrate 
the good numerical stability of this approach by computing six-, seven- and eight-gluon amplitudes in QCD and 
comparing against known analytic results. 



1. Introduction 

Quantitatively reliable predictions for back- 
ground processes will play an important role in 
ferreting out signals of new physics in experi- 
ments at the upcoming Large Hadron Collider 
(LHC) . New physics beyond the Standard Model 
is expected to emerge in these TeV-scale exper- 
iments. Known-physics backgrounds from elec- 
troweak, QCD, and mixed processes will also con- 
tribute events that may overwhelm or mimic new- 
physics signals. Uncovering and understanding 
the new-physics signals will require use of elabo- 
rate kinematic requirements (such as several iden- 
tified jets, cuts on missing transverse energy, etc.) 
and reliable knowledge of background processes. 

Leading-order calculations in QCD suffer from 
large uncertainties and therefore do not suffice 
to furnish the required quantitative knowledge of 
backgrounds. Next-to-leading order calculations 
are required jT]. Indeed, for a few signal, back- 
ground, or calibration processes (Higgs-boson 
production, top production, and distributions as- 
sociated with production of single electroweak 
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vector bosons), precision calculations at next-to- 
next-to-leading order (NNLO) are needed. For 
other processes, NLO will likely suffice until the 
LHC enters its precision-measurement era. There 
are a large number of processes that ought to be 
computed, however, and these include many pro- 
cesses with many final-state jets, corresponding 
to large final-state parton multiplicity. 

NLO predictions can be provided by parton- 
level Monte Carlo programs, or by parton show- 
ers such as MC@NLO [2 . Both types of pro- 
gram require the computation of virtual and real- 
emission amplitudes. The computation of the lat- 
ter relies on well- understood technologies [314] . 
Parton-level programs also require the isolation 
of infrared singularities in the integrated real- 
emission contributions, and a means of can- 
celing them systematically against the virtual- 
correction singularities. These technologies |5|6j 
are also well understood and recently authors 
have begun to automate them [7j. The infrared- 
divergent parts of the one-loop virtual correc- 
tions are also well-understood |5|8j . The remain- 
ing, infrared-finite parts of these one-loop am- 
plitudes have been the difficult part of compu- 
tations with three final-state objects, and have 
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been the primary bottleneck to computations of 
processes with four or more final-state objects. 
BlackHat is one of a new generation of numeri- 
cal approaches that aims to break this bottleneck. 
Other numerical efforts along similar lines are de- 
scribed in refs. [9110111] . 

The traditional approach to one-loop computa- 
tions uses Feynman diagrams. With an increasing 
number of external particles, however, the com- 
putational complexity of such computations with 
traditional methods grows factorially. This com- 
plexity cannot be tamed by use of technologies 
such as the spinor-hclicity formalism [12] alone. 
Recent years have witnessed a ferment of devel- 
opment, based on the analytic properties of uni- 
tarity and factorization that any amplitude must 
satisfy, of new approaches to overcoming these 
computational difficulties [13|14|15I16|T7TL8] . in- 
cluding on-shell methods [19|20|21|22|23|24|25[ 
I26l27l28l29l30l3ll32l33l34l35l36| . These methods 
are efficient, and feature only mild growth in re- 
quired computer time with increasing number of 
external particles. They effectively reduce loop 
calculations to tree-like ones, simplifying them 
intrinsically and further allowing use of efficient 
algorithms for the tree- amplitude ingredients. 

2. The On-Shell Approach 

BlackHat is built using an on-shell ap- 
proach, the unitarity method [19138] with multi- 
ple cuts [5D] (a.k.a. generalized unitarity), along 
with significant refinements [22 30 33 exploiting 
complex momenta. The cuts replace two or more 
propagators by delta functions, thus putting the 
corresponding momenta on shell. They thereby 
isolate terms in the amplitude with distinct an- 
alytic structure, allowing them to be computed 
independently. 

We employ a predominantly numerical ap- 
proach in BlackHat. The loop integrals are 
evaluated from analytic formulae, likewise their 
contributions to residues at spurious singulari- 
ties (see section [4]); and analytic formulae may be 
used to speed up evaluation of tree amplitudes. 
Otherwise, the code is numerical, and in partic- 
ular everything that corresponds to algebra in a 
symbolic calculation is done numerically. This 



makes it easier to design a general-purpose code, 
as distinct from the bespoke analytic process-by- 
process calculations that have been done to date. 
It is also important, however, to obtain an effi- 
cient code. As the example of the Berends-Giele 
recursion relations shows [37] , evaluating them re- 
cursively with caching of intermediate currents, a 
numerical approach can make it much more prac- 
tical to eliminate repeated evaluation of common 
subexpressions than an analytic approach. In- 
deed, it seems likely that only a numerical ap- 
proach can meet the goal of a polynomial-time 
algorithm for the evaluation of each helicity am- 
plitude at one loop. 

We begin by separating the amplitude into cut- 
containing and rational parts, 



-<*-n — C n ~r -fin i 



(1) 



where the former contain all (poly) logarithms, 
■k terms, and the finite constant in the scalar 
bubble. The cut-containing part of massless 
dimensionally-regulated amplitudes with four- 
dimensional external momenta may be written in 
a basis of scalar integrals |39I40I41I42|43] . 



L ,1 



Y J dj\+Y. c ^ + Y. b * 1 ^ 



(2) 



The integrals -(2,3,4 are respectively bubble, trian- 
gle, and box integrals. 

We compute the coefficients bi, Ci, and di nu- 
merically using the unitarity method, with four- 
dimensional loop momenta. The tree amplitudes 
that feed into the computation may be evaluated 
efficiently using spinorial methods. We compute 
the rational parts using loop-level on-shell recur- 



3. Cut Parts 

The most straightforward coefficients to obtain 
are those of the box integrals. Imposing four cut 
conditions on a one-loop integrand, and then set- 
ting e = 0, freezes it completely. Furthermore, 
this isolates the coefficient of a single box inte- 
gral uniquely. The coefficient is then given [22j in 
terms of the product of four tree amplitudes at 
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the corners of the box, 



ja a tree a tree /(tree a tree 
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where the cut loop momenta /,■ are the two so- 
lutions to the quadruple-cut on-shell conditions, 
labeled by o = ±. 

In the case of triangle coefficients, we can im- 
pose three cut conditions. Here, the cut condi- 
tions no longer freeze the integrand completely; 
one degree of freedom is left over. There are dif- 
ferent ways to parametrize this degree of free- 
dom [14130] ; we use the variant proposed by 
Forde [331, 



t 
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where K i t3 are the two (possibly massive) exter- 
nal momenta separated by the cut propagator I; 
1 = 1± = -K 1 -K 3 ± VA; Si = Kf, and 
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Once the remaining degree of freedom is 
parametrized by t, the integrand has the following 
form, 
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The sum over poles corresponds to the vari- 
ous box contributions sharing the same triple 
cuts. We follow Ossola, Papadopoulos, and Pit- 
tau (OPP) [3D] in subtracting these contributions 
from the integrand, leaving behind seven indepen- 
dent coefficients. Forde's parametrization isolates 
the desired coefficient cq by virtue of its analytic 



property — it is the constant as t — > 00, or equiv- 
alently it can be extracted as the residue of the 
integrand divided by t, 



Co 
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We evaluate this contour integral using a discrete 
Fourier projection, 
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where £0 is an arbitrary complex number. The 
projection avoids the potentially-unstable matrix 
inversion that could arise from simply inverting 
a system of equations to solve for the seven co- 
efficients c_3, . . . , c 3 . (The other coefficients are 
needed in order to subtract triangle coefficients 
when computing in their turn the bubble ones, 
and can also be obtained by a Fourier projection.) 
The subtraction of the box coefficients makes the 
discrete Fourier projection exact and also allows 
for the flexibility in the choice of to- It thereby 
improves the numerical stability of the calcula- 
tion. 

The bubble coefficients can be computed fol- 
lowing the same approach, subtracting both box 
and triangle contributions, and using a two- 
dimensional discrete projection. We refer the 
reader to ref. [H] for more details. 

4. Rational Parts 

On-shell recursion relations for the rational 
terms may be derived by considering deforma- 
tions of the amplitude, parametrized by a com- 
plex parameter z 24J. These deformations shift 
two external momenta by ±z • q where q is a com- 
plex null four- vector, so as to preserve overall mo- 
mentum conservation and leave all external mo- 
menta on shell. The recursion relation follows 
from evaluating the contour integral 
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where C is a circle at infinity. If i?Jf rgc z does 
not vanish for a given choice of deformation, it 
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Figure 1. The distribution of the logarithm of the relative error over 100,000 phase- 

space points for the MHV amplitudes A 6 (l~ , 2",3 + ,4+, 5+, 6+), A 7 (l~ , 2~,3+,4+, 5+, 6+, 7+) and 
A 8 (l~,2~,3 + ,4 + ,5 + ,6 + , 7 + ,8 + ). The dashed (black) curve in each histogram gives the relative er- 
ror for the 1/e 2 part, the solid (red) curve gives the 1/e singularity, and the shaded (blue) distribution 
gives the finite e° component of the corresponding helicity amplitude. 



can be computed using an auxiliary recursion re- 
lation [27] . The rational terms may be computed 
using Cauchy's theorem, 



R n (0) = R 



larti 



E 

poles a 



Res 



R n {z) 
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The sum over the poles in the last term can be de- 
composed into two sets, the physical or spurious 
poles, depending on whether the pole is or is not 
present in the full deformed one-loop amplitude, 



±\n Jl,, 
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The contributions R® from the physical poles 
can be computed using on-shell recursive dia- 
grams [26] . 

In BlackHat, the residues at the spurious 
poles are computed using the cut parts. Because 
the spurious poles must cancel in the amplitude 
as a whole, we have 



R n = - E 



Res 

z /3 
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where C n (z) is the cut part from eq. |T]), as de- 
formed by the on-shell deformation parametrized 
by z. 

The spurious singularities arise from zeros of 
Gram determinants implicitly appearing in the 
denominators of the integral coefficients bi, Cj, 
and di of eq. ^. We evaluate the residues of these 
singularities by making a discrete approximation 
to a contour integral on a small circle around the 
pole. At each complex value on the circle, we 
evaluate the coefficients bi, Ci, and di numerically 
as described in the previous section. Since the 
residues are of course rational, and can only arise 
from expanding the dilogarithms or logarithms 
in the integral functions, we do not directly eval- 
uate the loop integrals numerically; rather, we 
first expand them analytically in the appropriate 
neighborhood of the Gram-determinant singular- 
ity, and then evaluate the rational expansion co- 
efficients numerically. 

We can control the approximation by choosing 
the size of the circles around the spurious sin- 
gularity and the number of points on the circle. 
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(The current code evaluates the integrand at ten 
points around the circle.) This risks numerical 
instabilities if the circle becomes too small. 

In order to handle possible numerical insta- 
bilities, we test for them dynamically, that is 
event by event and spurious pole by spurious 
pole. We test to see that the coefficient of the 
non-logarithmic 1/e singularity is reproduced cor- 
rectly. It is sufficient to test the bubble coeffi- 
cients, because they produce this singular term. 
Box and triangle coefficients are tested indirectly, 
because they are subtracted in order to compute 
the bubble coefficients. We also test to see that 
the spurious singularities cancel in the sum over 
bubble coefficients. If either of these tests fail, 
we recompute the numerically unstable terms of 
the amplitude at higher precision. Because the 
fraction of points failing the tests is small, this 
does not impose a significant time penalty on our 
code. 



5. Results 

As an example, we used BlackHat to com- 
pute the one- loop 2 — > 4-, 2 — ► 5-, and 2^6- 
gluon maximally helicity-violating (MHV) ampli- 
tudes for nf = at 100,000 phase-space points, 
generated using a flat distribution. (We impose 
the following cuts: Et > 0.01y/s, n < 3, and 
Afl > 0.4.) We compared the numerical results 
against those computed using known analytic re- 
sults. The histogram in fig. Q] shows the results. 
The horizontal axis gives the logarithmic relative 
error, 



log- 



in 



^target 



\A\ 



(14) 



for each of the I /e 2 , 1/e, and e° parts of the one- 
loop amplitude. The vertical axis in these plots 
shows the number of phase-space points in a bin 
that agree with the target to a specified relative 
precision. The vertical scale is logarithmic, which 
enhances the visibility of the tail of the distribu- 
tion, and thereby illustrates the good numerical 
stability of the computation. 
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